using System;
using System.Collections.Generic;
using System.Linq;

using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Plotting.Functions;

namespace Atomic.Thermodynamics.Electrons
{
	public class ElectronicPotentialOLD : IPotential
	{
		private IElectronicStateDensity stateDensity;
		[NonSerialized] private InnerEvaluator evaluator;
		[NonSerialized] private RollingCache<double, Tuple<double, double, double>> values;

		public ElectronicPotentialOLD(IElectronicStateDensity stateDensity)
		{
			this.stateDensity = stateDensity;
		}

		public double FreeEnergy(double temperature)
		{
			// Compute for each temperature. No linear interpolation! Plotting is very easy now, e.g. t => p.FreeEnergy(t).
			return Compute(temperature).Item1;
		}

		public double Entropy(double temperature)
		{
			return Compute(temperature).Item2;
		}

		public double HeatCapacity(double temperature)
		{
			return Compute(temperature).Item3;
		}

		/// <summary>
		/// The interpolated density of states function used for the calculations.
		/// </summary>
		public double StateDensity(double energy)
		{
			IPlotFunction density = new LinearInterpolationFunction(stateDensity.Points.Select(p => new PlaneVector(p.Energy - stateDensity.FermiEnergy, p.Density)));
			return density.Value(energy - stateDensity.FermiEnergy);
		}

		private Tuple<double, double, double> Compute(double temperature)
		{
			if (values == null)
			{
				values = new RollingCache<double, Tuple<double, double, double>>(1000);
			}

			if (evaluator == null)
			{
				evaluator = new InnerEvaluator(stateDensity);
			}

			Tuple<double, double, double> value;
			if (!values.TryGetValue(temperature, out value))
			{
				double f, s, cv;
				evaluator.Compute(temperature, out f, out s, out cv);

				// Cache computed free energy, entropy, and heat capacity as a tuple.
				value = Tuple.Create<double, double, double>(f, s, cv);
				values.Add(temperature, value);
			}

			return value;
		}

		[Serializable]
		private class InnerEvaluator
		{
			private IPlotFunction density;
			private int n;
			private double e1, e2, de;
			private double[] energies, densities;

			public InnerEvaluator(IElectronicStateDensity stateDensity)
			{
				// Use a IPlotFunction to perform linear interpolation.
				density = new LinearInterpolationFunction(stateDensity.Points.Select(p => new PlaneVector(p.Energy - stateDensity.FermiEnergy, p.Density)));

				n = 2501*5;
				e1 = -20.21425944;
				e2 = 11.92309043;
				de = (e2 - e1) / n;

				energies = new double[n];
				densities = new double[n];
				for (int i = 0; i < n; i++)
				{
					double e = e1 + i * de;
					energies[i] = e;
					densities[i] = density.Value(e);
				}
			}

			public void Compute(double t, out double f, out double s, out double cv)
			{
				double mu = FindFermiLevel(t);
				double kb = 1.380658e-23 / 1.60217733e-19;  // k_B in eV/K;

				double int1 = 0.0;
				double int2 = 0.0;
				for (int i = 0; i < n; i++)
				{
					double e = energies[i];
					double ff = FermiDistribution(e, t, mu);
					double d = densities[i];
					int1 += e*ff*d*de;
					if (e < mu)
					{
						int2 += e*d*de;
					}
				}
				double resE = int1 - int2;

				// Electronic entropy.
				double maxE=mu + 20*kb*t;
				double minE=mu - 20*kb*t;
				double dE=(maxE-minE)/200;

				s = 0.0;
				cv = 0.0; // heat capacity
				for (int i = 0; i < 200; i++)
				{
					double e = minE + dE * i;
					double ff = FermiDistribution(e, t, mu);
					double intSe = density.Value(e) * (ff * Math.Log(ff) + (1.0 - ff) * Math.Log(1.0 - ff));
					s -= kb * intSe * dE;
					cv += density.Value(e) * Math.Exp((e - mu) / (kb * t)) * Math.Pow(ff * (e - mu) / t, 2.0) / kb * dE; // from T*dS/dT
				}

				// Electronic free energy.
				f = resE - t * s;

				// FIXME check!
				//f /= 48; // eV (per atom)
				//cv /= 48;
			}

			private double FermiDistribution(double e, double t, double mu)
			{
				double kb = 1.380658e-23 / 1.60217733e-19;  // k_B in eV/K;
				return 1.0 / (Math.Exp((e - mu) / (kb * t)) + 1.0);
			}
			
			private double FindFermiLevel(double t)
			{
				double n0 = FindDensityIntegral(0.0, 0.5);

				// Determine Fermi energy using bisection. The interval is from Shang's original implementation.
				return RootFinder.Bisection(x => FindDensityIntegral(x, t) - n0, -0.6, 0.9, 20);
			}

			private double FindDensityIntegral(double mu, double t)
			{
				//double kb = 1.380658e-23 / 1.60217733e-19;  // k_B in eV/K;

				double res = 0.0;
				double y0 = 0.0;
				for (int i = 0; i < n; i++)
				{
					double e = energies[i];
					double fermi = FermiDistribution(e, t, mu);
					double y = fermi * densities[i];
					res += 0.5 * (y + y0) * de;
					y0 = y;
					//res += y * delta;
				}

				return res;
			}
		}
	}
}
